An efficient algorithm for computing the 

Baker-Campbell-Hausdorff series and some of its 

O ■ applications 
O ■ 
(N 



o 

in 



(N 

d 



oo 

o 



Fernando Casas^* Ander Murua^^ 



Abstract 



We provide a new algorithm for generating the Baker-Campbell-Haus- 
' (— l'' dorff (BCH) series Z = log(e'''"e^) in an arbitrary generahzed Hall basis 

O |. of the free Lie algebra C{X, Y) generated by X and Y. It is based on 

the close relationship of C{X,Y) with a Lie algebraic structure of labeled 
rooted trees. With this algorithm, the computation of the BCH series up 
to degree 20 (111013 independent elements in C{X,Y)) takes less than 
15 minutes on a personal computer and requires 1.5 GBytes of memory. 
We also address the issue of the convergence of the series, providing an 
optimal convergence domain when X and Y are real or complex matrices. 
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1 Introduction 



^ ■ The Baker-Campbell-Hausdorff formula deals with the expansion of Z in e^ e^ = 

^ ! in terms of nested commutators of X and Y when they are assumed to be 

non-commuting operators. If we introduce the formal series for the exponential 
function 

oo 

e^e^ = y —XPY" (LI) 

p,q=0^ ^ 

and substitute this series in the formal series defining the logarithm function 



iogz=j:i^(z-i)^- 



k 

k=l 



one obtains 



Iog(e e ) = 2^^— 2^ 
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where the inner summation extends over all non-negative integers pi, qi, . . . , 
Pk, Qk for which Pi + qi > {i = 1, 2, . . . , fc). Gathering together the terms for 
which Pi + (7i + P2 + 92 + ■ ■ ■ + Pfc + ^fc = we can write 



Z = log(e^e^) = (1.2) 



m=l 

where Pm{X, Y) is a homogeneous polynomial of degree m in the non-commuting 
variables X and Y . Campbell [8], Baker [2] and Hausdorff [17j addressed the 
question whether Z can be represented as a series of nested commutators of X 
and Y , without producing a general formula. We recall here that the commu- 
tator [X, Y\ is defined as XY — YX. It was Dynkin [12j who finally derived an 
explicit formula for Z as 

z-yy t^'"'^'' ■ • • ^'"^'"^ (13) 

Here the inner summation is taken over all non-negative integers pi,qi, . . pk, 
qk such that pi + qi > . . . ,pk + qt > ^ and [X^iyi . . . X^^y^''] denotes the 
right nested commutator based on the word X'^^ Y'^^ . . . X^^ Y'"' . Expression 
(jl.Sp is known, for obvious reasons, as the Baker-Campbell-HausdorjJ series in 
the Dynkin form. By rearranging terms, it is clear that Z can be written as 

oo 



log(e^ e^) = X + y + ^ Z^, (1.4) 



with Zjn{X, Y) a homogeneous Lie polynomial in X and Y of degree m, i.e., it is 
a Q-linear combination of commutators of the form [Vi, [V2) • • • , [Vm-i, Vm] ■ ■ ■]] 
with Vi £ {X, Y} for 1 < i < m. The first terms read explicitly 

Z2 = l[X,Y] 

Z, = 1[X,[X,Y]]-^[Y,[X,Y]] 
Za = ^[X,[Y,[Y,X]]] 

The expression e^ e^ = e^ is then called the Baker-Camphell-Hausdorff for- 
mula (BCH for short), although other different labels (e.g., Campbell-Baker- 
Hausdorff, Baker-Hausdorff, Campbell-Hausdorff ) are commonly attached to 
it in the literature. The formula (jl.Sp is certainly awkward to use due to the 
complexity of the sums involved. Notice, in particular, that different choices 
oi Pi, qi, k in (11. 3p may lead to terms in the same commutator. Thus, for 
instance, [X'^Y'^] = [X^'Y^X'^Y^] = [X,[X,[X,Y]]]. An additional difficulty 
arises from the fact that not all the commutators are independent, due to the 
Jacobi identity [17] : 

[Xi, [X2,X3]] + [X2, [Xs,Xi]] + [X3, [Xi,X2]] = 0. 
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The BCH formula plays a fundamental role in many fields of mathematics 
(theory of linear differential equations [26] , Lie groups [H] , numerical analysis 
[IB]), theoretical physics (perturbation theory [TU], Quantum Mechanics [15] . 
Statistical Mechanics [Ml [SO], quantum computing [ID]) and control theory 
(analysis and design of nonlinear control laws, nonlinear filters, stabilization of 
rigid bodies [IS]). In particular, in the theory of Lie groups, with this formula 
one can explicitly write the operation of multiplication in a Lie group in canon- 
ical coordinates in terms of the Lie bracket operation in its tangent algebra and 
also prove the existence of a local Lie group with a given Lie algebra |14j . 

Also in the numerical treatment of differential equations on manifolds \19\ 
I16j . the BCH formula is quite useful. If is a smooth manifold and X{M) 
denotes the linear space of smooth vector fields on A4, then a Lie algebra 
structure is established in X{^A) by using the Lie bracket [X, Y] of fields X and 
Y G XiM) [17] . The fiow of a vector field X G X{M) is a mapping exp{X) 
defined through the solution of the differential equation 



as exp{tX)(q) = u(t). Many numerical methods used to approximately solv- 
ing equation (II. 5p are based on compositions of maps that are flows of vector 
fields [16]. To be more specific, suppose the vector field X can be split as 
X = A + B and that the flows corresponding to A{u) and B{u) can be ex- 
plicitly obtained. Then one may consider an approximation of the form = 
exp(/iaiA) ex.p{hbiB) • • • exp{hakA) exp{hbkB) for the exact flow exp{h{A + B)) 
of (jl.Sp after a time step h. The idea now is to obtain the conditions to be sat- 
isfied by the coefficients Oj, bi so that */i((?) = u{h) + 0{hP~^^) as /i ^ 0, and 
this can be done by applying the BCH formula in sequence to the expression 
of "if up to the degree required by the order of approximation p [27]. This 
task can be carried out quite easily provided one has explicit expressions of 
implemented in a symbolic algebra package [23l HB] . 

In addition to the Dynkin form ()1.3p . there are other standard procedures 
to construct explicitly the BCH series. Recall that the free Lie algebra C{X, Y) 
generated by the symbols X and Y can be considered as a subspace (the sub- 
space of Lie polynomials) of the vector space spanned by the words w in the 
symbols X and Y, i.e., w = aia2 ■ ■ ■ am, each ai being X or Y. Thus, the BCH 
series admits the explicit associative presentation 



in which g^j is a rational coefficient and the inner sum is taken over all words 
w with length \w\ = m. Here the length of w is the number of letters it 
contains. The coefficients can be computed with a procedure based on a family 
of recursively computable polynomials |13j . 

Although the terms in equation (jl.Op are expressed as linear combinations 
of individual words (which are not Lie polynomials) , by virtue of the Dynkin- 




u{0) = qeM 



(1-5) 



oo 




(L6) 
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Specht-Wever theorem |21j . Z can be written as 

oo 

Z = X + F+^ - E (1-7) 

m=2 w\w\=m, 

that is, the individual terms are the same as in the associative series p.6p except 
that the word w = 0102 .. . am is replaced with the right nested commutator 
[w] = [ai,[a2, . . . [flm-i) o,m] ■ ■ ■]] and the coefficient Qw is divided by the word 
length m [12]. This gives explicit expressions of the terms Zm in the BCH 
series (|1.4p as a linear combination of nested commutators of homogeneous 
degree, that is, as a linear combination of elements of the homogeneous subspace 
C{X,Y)m of degre m of the free Lie algebra C{X,Y). However, it should be 
stressed that the set of nested commutators [w] for words w of length m is not 
a basis of the homogeneous subspace C{X,Y)m- 

By introducing a parameter r and differentiating with respect to r the power 
series '}2m>i ^"^^m = log(exp(rX) exp(Ty)), the following recursion formula is 
derived in [47] : 

Zi = X + Y (1.8) 

, [{-"l)/2] 

mZm = -[X-Y,Zm-i]+ Yl 7^^'^"^^ + ^))^' 

p=i 

Here Z = J2ni>i ^rn, ad|(X + = [Z,ad^-^(X + Y)], the stand for the 

Bernoulli numbers [1], and (ad^(X + 5^) ) denotes the projection of ad^(X + 

Y) onto the homogeneous subspace C{X,Y)fn, which can be written in terms 
of Zi, Z2, Z3, ... as 

Ud'iiX + Y)) = [Z,,,[---[Zk,,,X + Y]-..]]. 

fcj^H \'k2p—ffi'~^ 

fel>l,...,fc2p>l 

Explicit formulas (jl.3p . (jl.Tp as well as recursion (jl.Sp can be used in princi- 
ple to construct the BCH series up to arbitrary degree in terms of commutators. 
As a matter of fact, several systematic computations of the series have been car- 
ried out along the years, starting with the work of Richtmyer and Greenspan in 
1965 [37] ■ where results up to degree eight are reported. Later on, Newman and 
Thompson obtained the coefficients in (jl.7p up to words of length 20 [32] , 
Bose [6] constructed an algorithm to compute directly the coefficient of a given 
commutator in the Dynkin presentation (|1.3p and Oteo [33] and Kolsrud [22] 
presented a simplified expression of (11. 3p in terms of right nested commutators 
up to degree eight and nine, respectively. More recently, Reinsch [35] has pro- 
posed a matrix operation procedure for calculating the polynomials Pm.{X, Y) in 
(jl.2p which can be easily implemented in any symbolic algebra package. Again, 
the Dynkin-Specht-Wever has to be used to write the resulting expressions in 
terms of commutators. 

As mentioned before, all of these procedures exhibit a key limitation, how- 
ever: the iterated commutators are not all linearly independent due to the 
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Jacobi identity (and other identities involving nested commutators of higher 
degree which are originated by it [33] )• In other words, they do not provide 
expressions directly in terms of a basis of the free Lie algebra C{X,Y). This 
is required, for instance, in applications of the BCH formula in the numerical 
integration of ordinary differential equations, or when one wants to study spe- 
cific features of the series, such as the distribution of the coefficients and other 
combinatorial properties |32j . 

Of course, it is always possible to express the resulting formulas in terms 
of a basis of C{X, Y) but this rewriting process is very time consuming and 
requires a good deal of memory resources. In practice, going beyond degree 
m = 11 constitutes a difficult task indeed [281 [23| I46j . since the number of 
terms involved in the series grows exponentially. 

Our goal is then to express the BCH series as 

Z = log(exp(X) exp(y)) = Y,Zi E^, (1.9) 

where G Q (i > 1) and {Ei : i = 1, 2, 3, . . .} is a basis of C{X,Y) whose 
elements are of the form 

Ei = X, E2 = Y, and Ei = [Ei> , Ei>i] i>3, (1.10) 

for appropriate values of the integers < i (for i = 3, 4, . . .). Clearly, each 
Ei in (jl.lOp is a homogeneous Lie polynomial of degree where 

|1| = |2| = 1, and \i\ = + for i > 3. (1.11) 

We will focus on a general class of bases of the free Lie algebra C{X, Y), referred 
to in the current literature as generalized Hall bases and also as Hall-Viennot 
bases [361 SH] • These include the Lyndon basis [251 SH] , and different variants 
of the classical Hall basis (see [36] for references). Specifically, in this paper 
we present a new procedure to write the BCH series (ll.9p for an arbitrary 
Hall-Viennot basis. Such an algorithm is based on results obtained in |30j . in 
particular those relating a certain Lie algebra structure g on rooted trees with 
the description of a free Lie algebra in terms of a Hall basis. This Lie algebra g 
on rooted trees was first considered in [TT] , whereas a closely related Lie algebra 
on labeled rooted trees was treated in [15] (see [18] for the relation of these two 
Lie algebras and for further references about related algebraic structures on 
rooted trees). 

We have implemented the algorithm in Mathematica (it can also be pro- 
grammed in Fortran or C for more efficiency). The resulting procedure gives 
the BCH series up to a prescribed degree directly in terms of a Hall-Viennot 
basis of C{X,Y). As an illustration, obtaining the series (in the classical ba- 
sis of P. Hall) up to degree m = 20 with a personal computer (2.4 GHz Intel 
Core 2 Duo processor with 2 GBytes of RAM) requires less than 15 minutes 
of CPU time and 1.5 GBytes of memory. The resulting expression has 109697 
non-vanishing coefficients out of 111013 elements Ei of degree \i\ < 20 in the 
Hall basis. As far as we know, there are no results up to such a high degree 
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reported in the literature. For comparison with other procedures, the authors 
of [l6] report 25 hours of CPU time and 17.5 MBytes with a Pentium III PC 
to achieve degree 10. By contrast, our algorithm is able to achieve m = 10 in 
0.058 seconds and only needs 5.4 MBytes of computer memory. 

In Table [3] in the Appendix, we give the values of i' and i" for the elements 
Ei of degree |i| < 9 in the Hall basis and their coefficients Zi in the BCH formula 
(jl.Op . The elements of the basis are ordered in such a way that i < j if |z| < 
and the horizontal lines in the table separate elements of different homogeneous 
degree. Extension of Table [3] up to terms of degree 20 is available at the website 
www.gicas.uji.es/research/bch.html for both the basis of P. Hall and the 
Lyndon basis. As an example, the last element of degree 20 in the Hall basis is 

^111013 = mY,X],Y],[Y,X]],[[[Y,X],X],[Y,X]]], 
[[[[Y,X],Y],[Y,X]],[[[[Y,X],Y],Y],Y]]], 

and the corresponding coefficient in ()1.9p reads 

19234697 
2111013 - ~ 140792940288 ' 
Another central issue addressed in this paper concerns the convergence prop- 
erties of the BCH series. Suppose we introduce a sub-multiplicative norm || ■ || 
such that 

<^||x||||y|| (1.12) 

for some fi > 0. Then it is not difficult to show that the series (jl.3p is absolutely 
convergent as long as \\X\\ + ||y || < (log 2)//i OUT]. As a matter of fact, several 
improved bounds have been obtained for the different presentations. Thus, in 
particular, the Lie presentation (II. 7p converges absolutely if \\X\\ < and 
11^ II < I/a* in a normed Lie algebra g with a norm satisfying (jl.l2p |3H I45j . 
whereas in [3] it has been shown that the series Z = J2m>i is absolutely 
convergent for all X, Y such that 

mII^II < / n TTTrdt (1.13) 

y^iiyii 2 + 1(1 -COt(it)) 

and the corresponding expression obtained by interchanging in (jl.l3p X by 
Y. Moreover, the series diverges in general if ||X|| + ||y|| > vr when /_f = 2 
|28j . Here we provide a generalization of this feature based on the well known 
Magnus expansion for linear differential equations [26] and also we give a more 
precise characterization of the convergence domain of the series when X and Y 
are (real or complex) matrices. 

2 An algorithm for computing the BCH series based 
on rooted trees 

2.1 Summary of the procedure 

Our starting point is the vector space g of maps a : T ^ M, where T denotes 
the set of rooted trees with black and white vertices 
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In the combinatorial literature, T is typically referred to as the set of labeled 
rooted trees with two labels, 'black' and 'white'. Hereafter, we refer to the 
elements of T as bicoloured rooted trees. 

The vector space g is endowed with a Lie algebra structure by defining the 
Lie bracket [a, /3] G 3 of two arbitrary maps a, /3 S g as follows. For each u T, 

{a,P]{u) = - a(7x(^'))/3(n(,.))), (2.1) 

where |n| denotes the number vertices of u, and each of the pairs of trees 
{u(^j^,u^^^) € T X T, j = 1, . . . , — 1, is obtained from u by removing one of 
the \u\ — 1 edges of the rooted tree u, the root of being the original root of 
u. For instance, 

[a,/3](5) = a(o)/3(.)-«(.)/3(o), [a,/?](§) = 0, 

[«,/?]( V) = 2(a( J )/?(.)- «(.)/?( 5)), (2.2) 
[a,/?](V) = «(§)/?(. ) + «(5)/3(o)-a(.)/3(§)-«(o)/3(5). 

An important feature of the Lie algebra g is that the Lie subalgebra of g 
generated by the maps X,Y G Q defined as 

1 if n = « \ f 1 if ti = o 



is a free Lie algebra over the set {X, y} [30]. In what follows, we denote as 
C{X, Y) the Lie subalgebra of g generated by the maps X and Y. 

It has also been shown in [30] that for each particular Hall-Viennot basis 
{Ei : i = 1, 2, 3, . . .}, (whose elements are given by (jl.lOp for appropriate 
values of i' , i" < i, i = 3,4, . . ., and X and Y given by (12. 3p ) one can associate a 
bicoloured rooted tree Ui to each element Ei such that, for any map a G /^{X, Y), 



i>i ^ ^' 



where for each i, a{ui) is certain positive integer associated to the bicoloured 
rooted tree Ui (the number of symmetries of Ui, that we call symmetry number 
ofui). For instance, the bicoloured rooted trees Uj and their symmetry numbers 
a{ui) associated to the elements Ei (of degree |i| < 5) of the Hall basis used in 
this work are displayed in Table [TJ 

As in the introduction, we denote by C{X,Y)n (n > 1) the homogeneous 
subspace of C{X, Y) of degree n (whence admiting {Ei : \i\ = n} as a basis). It 
can be seen [30] that, if a € C{X, Y), then its projection a„ to the homogeneous 
subspace C{X,Y)n is given by 

a (u) = i. "^^^ (2 5) 

1^ otherwise 

for each u (z T. 
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We also use the notation C{X, Y) for the Lie algebra of Lie series, that is, 
series of the form 



a = ai + Q2 + as H , where q„ e C{X, Y)n. 



Notice that in this setting, a Lie series a G C-{X, y) is a map a : T ^ R satisfy- 
ing that, for each n > 1, the map an given by (j2.5p belongs to C{X, Y)n- A map 
a € is then a Lie series if and only if (j2.4p holds (see [30] for an alternative 
characterization of maps a : T ^ M that actually belong to C{X,Y)). 

In particular, the BCH series Z = Zi + Z2 + + ■ ■ ■ given by (jl.Sp (for X 
and Y defined as in (|2.3|) ) is a Lie series. From (|1.8|) . it follows that Z( , ) = 
Z( o ) = 1, and for ?i = 2, 3, 4, . . . 



[{n-l)/2] 

(2p)! 



1'- ~ /„^2p. 



nZ(n) = -[X-Y,Z]{u)+ ^ -^(adf(X + y))(n) (2.6) 



p=i 

for each u (z T with n = Recall that, for arbitrary a,/3 G g and u the 
value [a, /3](m) is defined in terms of bicoloured rooted trees U(^j^,u^^^ with less 
vertices than n, so that (12. 6p effectively allows us to compute the values Z{u) 
for all bicoloured rooted trees with arbitrarily high number \u\ of vertices. In 
this way, the characterization (12. 4p of maps a G g that are Lie series directly 
gives a way to write Z S C{X, Y) in the form (jl.Op with 

z, = ^ for ^>1. (2.7) 

For instance, we have according to Table [U that in the Hall basis, 

i>l i>l '^^^'> 

= Z{.)X + Z{o)Y + Z{l)[Y,X] 

+^^^[%XlX]+Z{\f)[[Y,XlY] + -.. , 

where the first five coefficients Z[ui) can be obtained by applying (12. 6p with 
([22]): 

[x-y,z](j) = -z(.)-^(o) = -2, 

2Z(5) = i[x-y,z](5) = -1, 

[x-y,z](g) = 0, 2Z(§) = 0, 

[x-y,z](V) = -2^(5) = 1, 

[x-y,z](V) = ^(S)-^(§) = (2.8) 

[z,[z,x + y]]( V) = -2Z(.)(^(.)-^(o)) = 0, 

[z,[z,x + y]]( V) = -^(.)(^(.)-^(o)) = 0, 

3Z(V) = \^ 
3^(V) = "i- 
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In summary, the idea of the formahsm is to construct algorithmically a sequence 
of labeled rooted trees in a one-to-one correspondence with a Hall basis, ver- 
ifying in addition (j2.4p . In this way it is quite straightforward to build and 
characterize Lie series, and in particular, the BCH series. 



2.2 Detailed treatment 

In this subsection we provide a detailed treatment of the main steps involved in 
the procedure previously sketched, first by analyzing the representation (j2.4p 
of Lie series for the classical Hall basis and then by considering Hall-Viennot 
bases. 

We start by providing an algorithm that constructs the table of values (i' , i") 
(for z > 3) in (|1.10p (together with \i\ for i > 1) that determines a classical Hall 
basis. The algorithm starts by setting 

l' = l, 1" = 0, 2' = 2, 2" = 0, |1| = 1, |2| = 1, 

and initializing the counter i as i = 3. Then, the values \i\ for subsequent 
values of i are set as follows (i+~'~ indicates that the value of the counter i is 
incremented by one): 

Algorithm 1 

for n = 2,3,... 

j = l,...,i-l 

k = j + 

If \j\ + \k\ = n and j > k" then 
i" = J, i' = k, \i\ = n, 

The values of z', i" , \i\ thus determined satisfy that i' > i" > {i')" for i > 3. In 
addition, j < i if |j| < which implies that i', i" < i for all i > 3. The values 
for i' , and i" and the element Ei of the basis for the values of the index i of 
degree |«| < 5 are displayed in Table [TJ 

On the other hand, it is possible to design a simple recursive procedure to 
define the bicoloured rooted trees Ui appearing in (j2.4p in terms of the values 
of i' and i" by using the following binary operation. Given u,v ^ T, the new 
rooted tree n o -y s T is a rooted tree with |u| + vertices obtained by grafting 
the rooted tree v to the root of u (that is to say, n o is a new bicoloured rooted 
tree with the coloured vertices of u and v, one edge that makes the root of f a 
child of the root of u added to the edges of u and v). For instance, 

S = , and also o , = . 

We now define 

ui = « , U2 = o , and Uj = Uj/ oUj// for i >3. (2-9) 
Finally, the symmetry numbers = a{ui) can also be determined recursively: 
o"! = (T2 = 1, and (Tj = Ki (Tiiaiii, for i > 3, (2.10) 
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1 


1 

360 



Table 1: First elements Ei of the basis of P. Hall, their corresponding bicoloured 
rooted trees Ui, the values |i|, i" , i\ a{ui), and the coefficients Zi = Z{ui) / a{ui) 
in the BCH series ^M) 

where Ki = 1 if {i')" 7^ i" , and Ki = k^/ + 1 if {i')" = i" . 

The bicoloured rooted trees Uj, their symmetry numbers cr(Mj), and the 
coefhcients Zi = Z{ui) / a{ui) in the BCH series (jl.Op are displayed in Tabled] 
for the first values of the index z, whereas in Table [3] given in the Appendix, the 
terms of the BCH series ()1.9p up to terms of degree 9 are given in compact form 
for the classical Hall basis by displaying the values of i', i" , and Zi = Z{ui) /u[ui) 
for each index i. 

This procedure can be extended indeed to Hall-Viennot bases. A set {Ei : 
i = 1,2,3,...} C C{X,Y) recursively defined as (jl.lOp with some positive 
integers i' ,i" < i {i = 3,4,...) is a Hall-Viennot basis if there exists a total 
order relation >- in the set of indices {1, 2, 3, . . .} such that i >- i" for alH > 3, 
and the map 

d:{3,4,...} {{j,k)eZ+ xZ+ : jyktj"}, (2.11) 

d{i) = {i',i") (2.12) 

(with the convention 1" = 2" = 0) is bijective. 

In [l8l [36] , Hall-Viennot bases are indexed by a subset of words (a Hall set 
of words) on the alphabet {x,y}. Such Hall set of words {wi : i > 1} can be 
obtained by defining recursively Wi as the concatenation Wi'Wi" of the words 
Wi' and Will, with wi = x and W2 = y- For instance, the Hall set of words Wi 
associated to the indices i = 1, 2, . . . , 14 in Tableware x, y, yx, yxx, yxy, yxxx, 
yxxy, yxyy, yxxxx, yxxxy, yxxyy, yxyyy, yxxyx, yxyyx. 
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For the classical Hall basis we have considered before, the map ()2.1ip has 
been constructed in such a way that the total order relation >- is the natural 
order relation in Z"*", i.e., > (notice that in [7] the total order is chosen as <). 

This is not possible, however, for the Lyndon basis. The Lyndon basis can 
be constructed as a Hall-Viennot basis by considering the order relation >~ as 
follows: i y j if, in lexicographical order (i.e., the order used when ordering 
words in the dictionary) , the Hall word Wi associated to i comes before than the 
Hall word wj associated to j. The Hall set of words {wi : i > 1} corresponding 
to the Lyndon basis is the set of Lyndon words, which can be defined as the set 
of words w on the alphabet {x, y} satisfying that, for arbitrary decompositions 
of w as the concatenation w = uv oi two non-empty words u and v, the word 
w is smaller than v in lexicographical order [48\ [25] . 

Now, the representation (|2.4p of a map a € C{X, Y) (and in particular, the 
BCH series (|1.9|) with (|2.7p ) for any Hall-Viennot basis can be stated as follows. 

Theorem 2.1 Given a total order relation >- in 7L^ and a bijection \2.11\l 
satisfying that i >- i" for all i > 3, then any map a S C{X, Y) admits the 
representation ^2.4\ ) for the Hall basis \1.10\) and the bicoloured rooted trees Ui 
and their symmetry numbers (Ji = (T{ui) recursively defined as li2.9\) and ^2.10\) . 

Theorem 12.11 can be proven as a corollary of Theorem 3 and Remark 17 in 
|30| . Actually, in [30] it is shown that (j2.4p holds for a different set T = 
{ui,U2, Us, . . .} of bicoloured rooted trees associated to a Hall basis, for which 
a{ui) = 1 for all i. However, the set of Hall rooted trees we consider here (which 
is the set T* considered in Remark 17 in has some advantages from the 
computational point of view. 

In Table O we display the elements Ei of the Lyndon basis with degree 
|i| < 5, the corresponding Lyndon words vui, the bicoloured rooted trees Ui, the 
values i" , i' , o'{ui), and the coefficients Zi = Z{ui)/a{ui) in the BCH series 

(US]). 

2.3 Practical aspects in the implementation 

An important ingredient in the whole procedure is the practical implementation 
of the Lie bracket [a, /3] of two Lie series a, /3 € C{X, Y) C Q, which we address 
next. Let us consider for each u the sequence 

S{u) = {(^x(i), nW), . . . , (tx(|,|_i), n(H-i))} (2.13) 

of pairs of bicoloured rooted trees used to define the Lie bracket [a, [3] in (j2.ip . 
For instance, 

^(X) = {(V>o),(. ,<),(.,<)}. 

It can be seen that the sequences S{u) satisfy the following recursion. If n = 
vow, where v,uj £ T, then, let p = \v\ — 1, q = \w\ — 1, and 

S{v) = {(t;(i), ^;(^)), . . . , {v(^p),v^^^)}, S{w) = {(t/;(i), u;^), . . . , {w(,),w^^^)}, 
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Table 2: First elements Ei of the Lyndon basis, their corresponding Lyndon 
words Wi and bicoloured rooted trees Uj, the values |i|, i", z', a^Ui), and the 
coefficients Zj = Z{ui)/a{ui) in the BCH series (II. 9p 



then 

= (?;(i),'i;^^)ou;), . . . , (t;(p),?;(P)ou;), , ?;ou;(^) ), . . . , {w(^g^,Vow''''^)}. 

(2.14) 

From the point of view of implementation, it is important to observe that, if 
one wants to compute the BCH formula (jl.9p up to terms of degree |i| < n, 
there is no need to compute all values Z{u) for u (z T with \u\ < n, since only 
the values (12. 7p of Z{u) for the rooted trees Ui associated to the elements Ei of 
the basis are used in ()1.9p . However, when applying the recursion (12. 6p for the 
bicoloured rooted trees Ui, i = 1, 2, 3, . . . (that we will call Hall rooted trees), 
one finds that Z{u) needs to be computed for some additional bicoloured rooted 
trees u. Observe for instance that in (j2.8p . the value Z{u) for n = § was needed 
in order to get the value Z( *^ ) from the recursion ()2.6p . This is due to the 
fact that in (j2.ip . the bicoloured rooted trees uq) need not be Hall rooted trees 
when n is a Hall rooted tree (u^^^ is however necessarily a Hall rooted tree in 
that case). 

The minimal set %i of bicoloured rooted trees u for which Z{u) needs to be 
computed in order to get the values of Z[ui) for Hall rooted trees with |i| < n 
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by using recursion (j2.6p can be determined by requiring that 

{ui : \i\<n}cfnCT and S{fn) d % x 

It can be seen that the subset 7^ of bicoloured rooted trees can be alternatively 
defined as follows: We say that a bicoloured rooted tree ?; € T is covered by 
u G T if either v can be obtained from u by removing some of its vertices and 
edges, or u = t;. For instance, the bicoloured rooted trees covered by the tree 
uii in Table [2] are 

., o, I, 2> V> V> V*' X' sX- 

Then, it can be seen that 7^ is the set of bicoloured rooted trees covered by 
some of the trees in {uj : |i| < n}. 

As a summary of this treatment, we next describe the main steps of the 
algorithm that we use to compute the BCH series up to terms of a given degree 
for an arbitrary Hall-Viennot basis. Let mAr be sum of the dimensions of the 
homogeneous subspaces C{X,Y)n for 1 < n < iV, and let rfiN be the number 
of bicoloured rooted trees in T^v (so that m^v < '^at)- We proceed as follows 
for a given N: 

1. Determine the values i', i" for each i = 1, . . . ,m]y such that the Ei given 
by (jl.lOp are the elements of degree \i\ < N oi the required Hall-Viennot 
basis. Algorithm [1] can be used in the case of the basis of P. Hall. We 
use a similar (although slightly more complex) algorithm for the general 
case. 

2. Determine the bicoloured rooted trees u € T/v together with the |nj — 1 
pairs of bicoloured rooted trees in S{u) recursively obtained by (j2.14p . 
Actually, we associate each bicoloured rooted tree in T/v with a positive 
integer, such that T/v = {ui : i = 1, 2, . . . , ttiat} (and {uj : i = 
1,2,... ,m/v} is the set of Hall trees of degree \i\ < N). Each S{ui) is 
then represented as a list of |i| — 1 pairs of positive integers. 

3. Represent the truncated versions of Lie series a (truncated up to terms 
of degree A^) as a list of rfiN real values {ai, . . . , am^) corresponding 
to (a(ni), . . . , a(umjv))- The Lie bracket 7 = [a, (3] of two Lie series 
can be implemented as a way to obtain the list (71, . . . ,7mjv) from the 
lists («!, . . . , amjv) ^■^d (/3i, . . . , PniN) ™ terms of the pairs of integers 
representing S{ui) for each i = 1, . . . , ttt-tv- 

4. Represent the truncated versions of BCH series Z (truncated up to terms 
of degree A^) as a list of m/v rational values {Zi, . . . , Zf^j^) corresponding 
to {Z{ui), . . . , Z(umjv)), which can be obtained by initializing that list as 
(1, 1, 0, . . . , 0) and applying (j2.6p repeatedly for n = 2, . . . , A^. 

It is worth noticing that the number of trees in Tn is different for differ- 
ent Hall-Viennot bases. For instance, for the basis of P. Hall, T20 has 724018 
bicoloured rooted trees, while for the Lyndon basis the set T20 has 1952325 
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bicoloured rooted trees. Due to this fact, the amount of memory and CPU 
time required to compute with our algorithm the BCH formula up to a given 
degree for the Lyndon basis is considerably larger than for the basis of P. Hall. 
Moreover, the number of non-zero coefficients Zi in the BCH formula differs 
considerably in both bases. For instance, there are 109697 non-vanishing coef- 
ficients Zi (out of 111013 elements Ei of degree \i\ < 20) in the BCH formula 
for the basis of P. Hall, while for the Lyndon basis the number of non- vanishing 
coefficients Zi is 76760. 

3 Optimal convergence domain of the BCH series 

3.1 The BCH formula and the Magnus expansion 

One particularly simple way of obtaining a sharp bound on the convergence 
domain for the BCH series consists in relating it with the Magnus expansion 
for linear differential equations. For the sake of completeness, we summarize 
here the main features of this procedure. 

Suppose we have the non-autonomous linear differential equation 

^ = A{t)U, U{0) = I, (3.1) 

where U{t) and A(t) are operators acting on some Hilbert space Ti. (in particular, 
nx n real or complex matrices). Then the idea is to express the solution U{t) 
as the exponential of a certain operator f2(t), 

U{t) =expn{t). (3.2) 

By substituting (j3.2p into (j3.ip , one can derive the differential equation satisfied 
by the exponent 

oo „ 

^' = E#d^(^(*))' m = o. (3.3) 

fc=0 

By applying Picard's iteration on ()3.3p . one gets an infinite series for ^l{t), 

oo 

n{t) = J2 ^m{t), (3.4) 
m=l 

whose terms can be obtained recursively from 

Qiit) = [ A{ti) dti, 
Jo 

^mit) = {&d^(^s)A{s))mds m>2. (3.5) 

Equations (j3.2p and (j3.4p constitute the so-called Magnus expansion for the 
solution of (j3.ip . whereas the infinite series (j3.4p with (j3.5p is known as the 
Magnus series. 
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Since the 1960s |l9], the Magnus expansion has been successfully applied as 
a perturbative tool in numerous areas of physics and chemistry, from atomic and 
molecular physics to nuclear magnetic resonance and quantum electrodynamics 
(see [1] and [5] for a review and a list of references). Also, since the work 
by Iserles and N0rsett |2Uj . it has been used as a tool to construct practical 
algorithms for the numerical integration of equation (j3.1|) , while preserving the 
main qualitative properties of the exact solution. 

In general, the Magnus series does not converge unless A is small in a suitable 
sense, and several bounds to the actual radius of convergence have been obtained 
along the years. Recently, the following theorem has been proved [9]: 

Theorem 3.1 Let us consider the differential equation U' = A{t)U defined in 
a Hilhert space Tt, dim 7^ < oo, with U{0) = I, and let A{t) be a bounded linear 
operator on Tl. Then, the Magnus series Q{t) = X^^Li with il^ given by 

\3. 5|] converges in the interval t G [0, T) such that 

T 

\\A(s)\\ ds < TT 

and the sum Q{t) satisfies exp Q{t) = U{t). The statement also holds when 7i is 
infinite- dimensional if U is a normal operator (in particular, if U is unitary). 
Here || • || stands for the norm defined by the inner product on Ti. 

Moreover, it has been shown that the convergence domain of the Magnus 
series provided by this theorem is the best result one can get for a generic 
bounded operator A[t) in a Hilbert space, in the sense that it is possible to find 
specific A{t) where the series diverges for any time t such that ||^(s)||(is > tt 
[Ml [9]. 

Now, given two operators X and Y , let us consider equation p.ip with 

Clearly, the exact solution at i = 2 is given by U{2) = e"''"e^. On the other 
hand, if we apply recurrence (|3.5p to compute U{2) with the Magnus expansion, 
U{2) = e"(2)^ we get ni{2) = X + Y and more generally 17„(2) = Z„ in ([0]) . 
In other words, the BCH series can be considered as the Magnus expansion 
corresponding to the differential equation ()3.ip with A{t) given by ()3.6p at 
t = 2. 

Since ||A(s)||(is = ||X|| + ||y||, Theorem l3.1l leads to the following bound 
on the convergence of the BCH series. 

Theorem 3.2 Let X and Y be two bounded elements in a Hilbert space Ti with 
diniTi > 2. Then the BCH formula in the form i.e., expressed as a series 

of homogeneous Lie polynomials in X and Y , converges when \\X\\ + ||y|| < vr. 

Of course, this result can be generalized to any set Xi, X2, . . . , Xk of bounded 
operators: the corresponding BCH series is convergent if \\Xi || + • • • + || < tt 
in the 2-norm. 
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Let us illustrate the result provided by Theorem 13.21 with a simple example 
involving 2x2 matrices. 



x = { : M. y = (' i I (3.7) 



Example 1. Given 

^ -a J ' V 

with a, /3 € C, a simple calculation shows that 

iog(e"<='') = ^ + T:^>'. 

which is an analytic function for |q;| < tt with first singularities at a = ztivr. 
Therefore, the BCH formula cannot converge if \a\ > vr, independently of /? 7^ 0. 
By taking the spectral norm, it is clear that = |a|, = so that the 
convergence domain given by Theorem 13.21 is \a\ + < vr. Notice that in the 
limit — > this domain is optimal. □ 

Generally speaking, however, the bound given by Theorem 13.21 is conserva- 
tive, i.e., the BCH series converges for larger values of ||X|| and \\Y\\. Thus, 
in the previous example, for any a and /? with |q| < vr and |a| + > vr, the 
BCH series also converges. One would like therefore to have a more realistic 
characterization of this feature. It turns out that this is indeed feasible for 
complex n X n matrices. 



3.2 Convergence for matrices 

3.2.1 Convergence determined by the eigenvalues 

For complex n x n matrices it is possible to use the theory of analytic matrix 
functions and more specifically, the logarithm of an analytic matrix function, 
in a similar way as in the Magnus expansion [9] , to characterize more precisely 
the convergence of the BCH series. 

To begin with, let us introduce a parameter e G C and consider the substi- 
tution {X,Y) I — > {eX,eY) into eq. (|1.1|) . It is clear that 

U{e) = e"^ e'^ 

is an analytic function of e, det 7^ and the matrix function Z{e) = 
log [/(e) is also analytic at e = 0. Equivalently, the series Z(e) is convergent 
for sufficiently small e. It turns out that the actual radius of convergence of 
this series is related with the existence of multiple eigenvalues of [/(e). Let us 
denote by pi(e), . . . , /3„(e) the eigenvalues of the matrix [/(e). Observe that 
[/(O) = /, so that pi(0) = • • • = Pn(0) = 1 and we can take the principal values 
of the logarithm, logpi(O) = ••• = log/9„(0) = 0. In essence, if the analytic 
matrix function [/(e) has an eigenvalue poi^o) of multiplicity I > 1 for a certain 
eo such that: (a) there is a curve in the e-plane joining e = with e = eo, and 
(b) the number of equal terms in log/?i(eo), log/92(£o)7 • • • Aogpii^o) such that 
Pki^o) = PO) k = 1, . . . , Hs less than the maximum dimension of the elementary 
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Jordan block corresponding to po, then the radius of convergence of the series 
Z{e) = Ylik>i^^ verifying expZ(e) = U{e) is precisely r = |eo| [9]. 

More specifically, we find first the values of the parameter e for which the 
characteristic polynomial det({7(e) — pi) has multiple roots and write them in 
order of non-decreasing absolute value, 

^(1) ^(2) ^(3) /oo>. 

Next, we consider the circle |e| = l^o^^l in the complex e-plane and denote by 
/>Q^^ an eigenvalue of [/(eq^"*) with multiplicity li > 1. Let e move along some 
fixed curve L from e = to e = e^^^ in the circle \e\ < Then it is clear 

that li eigenvalues Pj{e) will tend to p^^^ at e = Sq^^ If these points lie at 
e = on the same sheet of the Riemann surface of the function logz, and 
this is true for all (possible) multiple eigenvalues of Yt{e) at e = eQ^\ then 
is called a extraneous root. Otherwise, Sq^^ is called a non- extraneous root. 

By the analysis carried out in [51], when |e| < l^o^^l the numbers log pj{e) 
are uniquely determined as eigenvalues of the matrix Z{e) and this series is 
convergent. This is also true at |e| = |eQ^^ | if eg^^ is an extraneous root, since then 
the eigenvalues of Z{£) retain their identity throughout the collision process, so 
that we proceed to the next value in the sequence (|3.8p until a non-extraneous 
root is obtained. 

(2) 

Assume, for simplicity, that Eq is the first non-extraneous root, for which 
there exists an eigenvalue po of U{e) with multiplicity I > 1. Associated with 
this multiple eigenvalue po there is a pair of integers (p, q) defined as follows. 

The integer p is the greatest number of equal terms in the set of numbers 
logpi(eo), logp2(eo), • • • ,logp;(eo) such that pk{eo) = po, k = I, . . . ,1. 

The integer q is the maximum degree of the elementary divisors {p — po)^ 
of U{£()), i.e., the maximum dimension of the elementary Jordan block corre- 
sponding to Po- 

Under these conditions, it has been proved that, if p < q for the eigenvalue 
Po, then the radius of convergence of the series Z{e) = J2k>i^^ is precisely 
r = \eo\ [M]. 

Although in some cases withp > q the series Z{e) may converge at |e| = |eo| 
and the radius of convergence r is greater than |eo| (for instance, when X and 
Y are diagonal), this situation is exceptional in a topological sense, as explained 
in [SU pp. 65-66]. 



3.2.2 Examples 

In order to illustrate this result we next consider a pair of examples involving 
also 2x2 matrices. 

Example 1. The first example involves again the matrices X and Y given by 
(|3.7p . In this case 
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The first values of e for which there are multiple eigenvalues of U{e) are 

TT 

e = 0, £ = ±i—. 

a 

The first value, e = 0, is clearly an extraneous root, whereas the eigenvalues 
of the matrix ?7(e) move along the unit circle, one clockwise and the other 
counterclockwise from 

Pi,2(0) = 1 to pi,2(«7r/a) = -1 

when e varies along the imaginary axis from e = to e = in/a (the same 
considerations apply to the case e = —in /a). Then, obviously, p = 1 and 
g = 2, so that the radius of convergence of the series Z{e) is 

vr 

By fixing e = 1, we get the actual domain of convergence of the BCH series as 
|a| = vr, i.e., the same result as in section \37\] □ 

Example 2. Consider now the matrices 

1 j' ^=(o 
and X = aA, Y = aB, with a > 0. Then 

U{e) =( ^ , 2] (3.9) 

has multiple eigenvalues when = 0, Eq'^^ = iti^. As e varies along the 
imaginary axis from e = to e = , the eigenvalues of the matrix U [e] , 

Pi,2(e) = \^—e^± 

move along the unit circle, one clockwise and the other counterclockwise from 

pi,2(0) = 1 to px^2{ef) = -1 

Thus, pi(eo ) and P2{£q ) he on different sheets of the Riemann surface of the 

(2) 

function logz and therefore Eg is a non-extraneous root, with p = 1. Since 

(2) 

U{eQ ) 7^ — /, we have g' = 2, so that the radius of convergence of the series 
Z{e) is precisely 

r=|4^^| = -. (3.10) 

This result should be compared with the bound provided by the Magnus ex- 
pansion. Since \\A\\ = \\B\\ = 1, Theorem 13.21 guarantees the convergence of 
the BCH series in this case whenever 2a|e| < vr, or |e| < which, in view of 
(j3.10p . is clearly a conservative estimate. 
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We can also check numerically the rate of convergence of the BCH series in 
this example as a function of the parameter e. Let us denote by Z^^l the sum 
of the first N terms of the series, i.e., 

N 

zy^\e) = Y,Zn{e) 

n=l 

and compute, for a = 2 and different values of e, the matrix 

= C/(e)e-^'"'(^)-/, 

where U{e) is given by (j3.9p . If e belongs to the convergence domain of the 
BCH series for the matrices X and Y (i.e., |e| < 1), then Er{e) — > as iV ^ oo. 

First we take e = |. With A'^ = 10, the elements of are of order 10"'^, 
whereas adding five additional terms in the series, N = 15, the elements of E^. 
are approximately 10^^'^. 

Next we choose e = 0.9, i.e., a value near the boundary of the convergence 
domain. In this case with = 15 the convergence of the series does not 
manifest at all. In fact, a much larger number of terms is required to achieve 
significant results. Thus, for the elements of Ej. to be of order 10~® we need to 
compute = 150 terms of the BCH series, whereas with = 200 the elements 
of Ej. are of order 10"^''. The computations have been carried out with the 
recurrence (jl.Sp . □ 

As this example clearly shows, it is not always possible to determine ac- 
curately the convergence domain of the BCH series by computing successive 
approximations, since the rate of convergence can be slow indeed near the 
boundary. For this reason it could be of interest to design a procedure to apply 
in practice the characterization of the convergence in terms of the eigenvalues 
of the matrix U (e) analyzed in subsection 13.2.11 for matrices. 

This procedure could be as follows. Given two matrices X, y, take the 
product of exponentials 

U{e) = e'^ e"^ 

with e = re*^. Next, define a grid in the e-plane, for instance in polar co- 
ordinates {r,0), by Ar = rf/{n + 1), A6 = 27r/(m + 1) for two integers 
71,171 > 1 and a sufficiently large value rf > 1. Then, for each point in the 
grid, (rfc = kAr, 9i = IA9), = 1, . . . , n + 1, / = 0, 1, . . . , m, compute the corre- 
sponding matrix f7(e) and its eigenstructure, locating where there are multiple 
eigenvalues (within a prescribed tolerance). If some of these multiple eigenval- 
ues have a negative real part, there exists a point in the neighborhood where the 
conditions enumerated in subsection 13.2.11 are satisfied and therefore we have 
approximately located the value of e where the BCH series fails to converge. 
This approximation can be made more accurate by applying, for instance, New- 
ton's method. The actual radius of convergence will be given by the smallest 
number r found in this way. Finally, if r > 1, then obviously the BCH series 
corresponding to X and Y converges. 
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4 Some Applications 



As an illustration of the usefulness of the previous results, in this section we 
present two not so trivial applications of the formalism developed in section 2 
for constructing explicitly the BCH series up to arbitrarily high order. 

4.1 The symmetric BCH formula 

Sometimes it is necessary to compute the Lie series W defined by 

exp(^X) exp(y) exp(ix) = exp(H^). (4.1) 

This occurs, for instance, if one is interested in obtaining the order conditions 
satisfied by time-symmetric composition methods for the numerical integration 
of differential equations \52\ I39j . Two applications of the usual BCH formula 
gives then the expression of W in the Hall basis of C{X, Y). 

A more efficient procedure is obtained, however, by introducing a parameter 
r in (j4.ip such that 

W^(T)=l0g(e^^/2^^^rX/2) (4 2) 

and deriving the differential equation satisfied by W{t). From the derivative of 
the exponential map, one gets 

^ = ^ + E ^""^wX, WiO) = Y (4.3) 

n=2 

whence it is possible to construct explicitly W as the series W{t) = X^^q (''"); 
with 

Wi{t) = Xt + Y 

W2{t) = (4.4) 



Wiir) = r{adi^X)ids, l>3 

j=2 



where, in general, W2m = for m > 1. By following a similar approach as with 
equation (II. 8p in the usual BCH series in section 2, the recursion (14. 4p allows 
one to express W in (|4.ip as 

W = Y,WiEi. (4.5) 

The coefficients Wi of this series up to degree 9 in the classical Hall basis are 
collected in Table H] in the Appendix. As with the usual BCH series, the 
coefficients up to degree 19 in both Hall and Lyndon bases can be found at 
www . gicas . u j i . es/research/bch . html. 

With respect to the convergence of the series. Theorem (j3.2p guarantees 
that W is convergent at least when \\X\\ + < vr. 
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4.2 The BCH formula and a problem of R.C. Thompson 

In a series of papers j43l 1321 I45|. R.C. Thompson considered the problem of 
constructing a representation of the BCH formula as 

e^ = e^ with Z = SXS^^ + TYT~^ , (4.6) 

for certain functions 5 = S{X^ y) and T = T{X, Y) depending on X and Y. By 
using analytic techniques related with the Kashiwara-Vergne method, Rouviere 
|38j proved that a Lie series p{X, Y) exists such that 

and converges when X, Y are replaced by normed elements near 0, whereas the 
representation (|4.6p is global when both X and Y are skew-Hermitian matrices 

m- 

Thompson himself developed a computational technique for constructing 
explicitly the series p{X, Y) up to terms of degree ten. Although his results 
were not published, he pointed out that they furnished strong evidence of the 
convergence of the series p{X, Y) on the closed unit sphere in any norm for 

which ||[x,y]|| < ||X|| ||y|| 05]. 

With the aim of clarifying this issue and illustrating the techniques de- 
veloped in section 2, we proceed next to compute p{X,Y). Since p{X,Y) G 
C{X,Y), i.e., is a Lie series, it can be written as 

p{X,Y) = Y,PiE^, 

i>l 

where the elements Ei have been introduced in ()1.9p . and the goal is to deter- 
mine the coefficients pi. This can be accomplished as follows. From the well 
known formula e^Ve~^ = e^'^'^V, it is clear that 

Z = e^'^p(x,Y)X + e'"^p(-'^'-^'>Y. (4.8) 

Next we expand e^'^p'-^'^'> X and e'**^''(~"^'--^)y into infinite series as a linear com- 
bination of the Hall basis in C{X, Y) and match the resulting terms with the 
corresponding to the BCH series for Z. Then a recursive system of equations 
is obtained for the coefficients pi. 

It is in fact possible to get a closed expression for p{X, Y) up to terms Y^ 
by taking into account the corresponding formula of Z [33] . Specifically, from 

Z = X + —^^^Y modY\ (4.9) 

a simple calculation leads to 

p{X, Y) = f{adx)Y mod Y^ 
with the function f{z) given by 

1 1 15 1 2 143 o 1 4 

f(z) = h-e^ = z-\ z^-\ z^-\ z^-\ (4.10) 

■'^ ' l-e"" z 4 96 384 92160 122880 ^ ^ 
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Working in the classical Hall basis, the complete expression up to degree four 
reads 



p{X,Y) 



\y + —\Y,X] + —\\Y,XIX] + — \\Y, X] , Y] 




11 



[[[Y,X],Y],Y] 



23040 



i.e., the corresponding equations have a unique solution. This is not the case, 
however, at degree 5, where a free parameter appears, which can be chosen to 
be pio- Then 



As a matter of fact, if higher degrees are considered, more and more free pa- 
rameters appear in the corresponding solution. Thus, at degree 7 there are two 
additional parameters (for instance, p26 and pso), whereas at degree 8 pso and 
P52 can be chosen as free parameters. We conclude, therefore, that there are 
infinite solutions to the problem posed by Thompson depending on an increas- 
ing number of free parameters. An interesting issue would be to determine the 
value of these parameters in order to render the whole series convergent on a 
domain as large as possible. 

4.3 Distribution of coefficients in the Lyndon basis 

As we previously mentioned, there are noteworthy differences in the results 
obtained when the algorithm of section [2] is applied to the BCH series in the 
classical P. Hall basis and the Lyndon basis, particularly with respect to number 
of vanishing coefficients. In the basis of P. Hall there are 1316 zero coefficients 
out of 111013 up to degree m = 20, whereas in the Lyndon basis the num- 
ber of vanishing terms rises to 34253 (more than 30% of the total number of 
coefficients). 

More remarkably, one notices that the distribution of these vanishing coef- 
ficients in the Lyndon basis follows a very specific pattern. Before entering into 
the details, let us denote for simplicity Cm = C{X,Y)m- We first remark that, 
for each m > 2, the Lyndon basis Bm of Cm is a disjoint union Bm = Bm,i^Bm,2 
with Bm,2 = [X,Bm~i]- Thus, Cm. = Cm,i®Cm,2, where Cm,2 = [X,Cm-i], and 
I3m,k {k = 1,2) is a basis of Cm,k- In particular, ad^~^y € Bm- In this sense, 
from our computations we make two observations: First, the coefficient in the 
BCH formula of the element ad^~^y in the basis Bm is for even m. Second, 
the coefficients for the terms in Bm,i are also zero for even m. This gives a total 
number of 



vanishing coefficients of terms of degree m = 2p in the BCH formula written in 
the Lyndon basis. Thus, for instance, when p = 10, the number of total number 
of vanishing coefficients is nc(20) = dim(£2o) — dim(>Ci9) -|- 1 = 52377 — 27594-1- 
1 = 24784. 



Pl2 = 



137 - 184320^10 
184320 



Pn = 



511 - 737280pio 
737280 



nc(2p) = dim{C2p) - dim{C2p^i) + 1, 



p>2 
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With these considerations in mind, we can proceed next to explain the 
observed phenomena. First, notice that expression (j4.9p gives explicitly the 
last term of the BCH series in the Lyndon basis at each degree. By formally 
expanding in power series of adx we get 

oo 

Z = X + Y + -adxY + ad'^xY mod F^. 

k=2 

Since -B2n+i = for all n > 1, the coefficient of ad^y is non-vanishing only for 
odd values of k, or equivalently, for even values of the degree m. 

As for the remaining zero coefficients, let us consider at this point the sym- 
metric BCH formula ()4.ip again. Clearly the series (14. Sh only contains terms of 
odd degree, i.e., W = X^j>ol^2i+i5 where Wi £ Ci. By denoting P = X/2 and 
forming the composition exp(P) exp(Vl^) exp(— P) one gets trivially 

e^'e^e-^ = e^e^=e^, 

i.e., the standard BCH formula. In the terminology of dynamical systems, 
exp(Vl^) and exp(Z) are said to be conjugated. Alternatively, we can write 
exp(Z') = exp(adp) exp(VF), so that Z = exp(adp)l^. It is worth to write 
explicitly this relation for each term Zm G -Cm of the series Z = Ylm>o by 
separating the odd and even degree cases. Specifically, 

Z2p+1 = W2p+1 + ^ (2j)!)22i '^^^ ^2p-2i+l 
= E (2j - 1)!)22J^1 ^'^X 'W2P-2J+1. 

From these expressions, it is clear that Z2p+i contains terms in the whole sub- 
space £2p+i,i © J~'2p+i,2 (due to the presence of W2p+i), whereas Z2p belongs to 
the subspace C2p,2, whose dimension is equal to dim(£2p-i)- In other words, 
the remaining dim(£2p) — dim(/32p-i) must necessarily vanish. In this sense, 
the Lyndon basis seems the natural choice to get systematically the BCH series 
with the minimum number of terms. Nevertheless, compared to the basis of P. 
Hall, more CPU time and memory is required to compute the BCH with our 
algorithm in the Lyndon basis. In particular, 1.5 GBytes are required to com- 
pute the BCH formula up to degree 20 in the Hall basis, whereas 3.6 GBytes 
of memory are needed in the Lyndon basis. 

5 Concluding remarks 

The effective computation of the BCH series has a long history and is closely 
related with the more general problem of carrying out symbolic computations 
in free Lie algebras. In this work we have presented a new algorithm which 
allows us to get a closed expression of the series Z = log(e"^e^) up to degree 
20 in terms of an arbitrary Hall-Viennot basis of the free Lie algebra generated 
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by X and Y , C{X, y), requiring reasonable computational resources. As far as 
we know, no other results are available up to this degree in terms of a basis of 
C{X, Y). The algorithm is based on some more general results presented in [3D] 
on the connection of labeled rooted trees with an arbitrary Hall-Viennot basis 
of the free Lie algebra. 

We have carried out explicitly the computations to get the coefficients of the 
BCH series in terms of both the classical P. Hall basis and the Lyndon basis, 
with some noteworthy differences in the corresponding results, as analyzed in 
subsection 14.31 

We have also addressed the problem of the convergence of the series when 
X and Y are replaced by normed elements. In the particular case of X and 
Y being matrices, we have provided a characterization of the convergence in 
terms of the eigenvalues of e'^. 

Although here we have considered only the BCH series, it is clear that 
other more involved calculations can be done, as is illustrated for instance by 
the problem of R.C. Thompson studied in section W?2\ . As a matter of fact, we 
intend to develop a general purpose package to carry out symbolic computations 
in a free Lie algebra generated by more than two operators. 
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A Appendix 

In Table [3] we collect the indices i' and i' for i > 3 in (jl.lO|) for the classical 
Hall basis and the values of the coefficients Zi in the BCH formula (jl.Op up to 
degree 9, whereas in Table [H we gather the corresponding coefficients for the 
symmetric BCH formula ()4.ip . 
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Table 3: Table of values of i' and i' for z > 3 in (jl.lOp for the classical Hall 
basis and the values € Q in the BCH formula (11.91) . 



i 


r 


■It 

I 




i 


i! 


•II 

I 




i 


./ 

I 


•II 

V 


Zi 


1 






1 


44 


25 
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26 


2 
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3 
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16 
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5 
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94 
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9 


6 


1 
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52 


18 
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Table 4: Table of values of i' and i' for i > 3 in (jl.lOp for the classical Hall 
basis and the values Wi € Q in the symmetric BCH formula ()4.ip . 



i 


./ 

I 


■It 

I 


Wi 


i 


1! 


•II 

I 


Wi 


i 


■ / 

I 


■II 

I 


Wi 


1 


1 





1 


44 


25 


2 





87 


31 


3 


1/4608 


2 


2 





1 


45 


26 


2 





88 


32 


3 


23/134400 


3 


2 


1 





46 


27 


2 





89 


33 


3 


1/37800 


4 


3 


1 


-1/24 


47 


28 


2 





90 


34 


3 


1/23040 


5 


3 


2 


-1/12 


48 


29 


2 





91 


35 


3 


1/201600 


6 


4 


1 





49 


15 


3 





92 


15 


4 


193/6451200 


7 


4 


2 





50 


16 


3 





93 


16 


4 


53/483840 


8 


5 


2 





51 


17 


3 





94 


17 


4 


25/193536 


9 


6 


1 


7/5760 


52 


18 


3 





95 


18 


4 


1/22400 


10 


6 


2 


7/1440 


53 


19 


3 





96 


19 


4 


-13/1209600 


11 


7 


2 


1/180 


54 


20 


3 





97 


20 


4 


53/483840 


12 


8 


2 


1/720 


55 


21 


3 





98 


21 


4 


17/161280 


13 


4 


3 


1/480 


56 


22 


3 





99 


22 


4 


-3/44800 


14 


5 


3 


-1/360 


57 


9 


4 





100 


23 


4 


-19/322560 


15 


9 


1 





58 


10 


4 





101 


15 


5 


367/4838400 


16 


9 


2 





59 


11 


4 





102 


16 


5 


193/645120 


17 


10 


2 





60 


12 


4 





103 


17 


5 


247/604800 


18 


11 


2 





61 


13 


4 





104 


18 


5 


53/241920 


19 


12 


2 





62 


14 


4 





105 


19 


5 


1/33600 


20 


6 


3 





63 


9 


5 





106 


20 


5 


53/161280 


21 


7 


3 





64 


10 


5 





107 


21 


5 


193/403200 


22 


8 


3 





65 


11 


5 





108 


22 


5 


13/201600 


23 


5 


4 





66 


12 


5 





109 


23 


5 


-1/5600 


24 


15 


1 


-31/967680 


67 


13 


5 





110 


9 


6 


-1 1 / flf? A 't A A 

11/774144 


25 


15 


2 


-31/161280 


68 


14 


5 





111 


10 


6 


1/290304 


26 


16 


2 


-13/30240 


69 


7 


6 





112 


11 


6 


-1/15360 


27 


17 


2 


-53/120960 


70 


8 


6 





113 


12 


6 


-89/1209600 


28 


18 


2 


-1/5040 


71 


8 


7 





114 


13 


6 


-11/241920 


29 


19 


2 


-1/30240 


72 


42 


1 


127/154828800 


115 


14 


6 


-13/80640 


30 


9 


3 


-53/161280 


73 


42 


2 


127/19353600 


116 


9 


7 


1/12096 


31 


10 


3 


-11/12096 


74 


43 


2 


157/7257600 


117 


10 


7 


-1-1 1 A -\ c\ 

11/64512 


32 


11 


3 


-3/4480 


75 


44 


2 


367/9676800 


118 


11 


7 


1/33600 


33 


12 


3 


-1/10080 


76 


45 


2 


23/604800 


119 


12 


7 


-11/120960 


34 


13 


3 


-1/4032 


77 


46 


2 


79/3628800 


120 


13 


7 


1/35840 


35 


14 


3 


-1/6720 


78 


47 


2 


1/151200 


121 


14 


7 


-29/134400 


36 


6 


4 


-19/80640 


79 


48 


2 


1/1209600 


122 


9 


8 


211/1935360 


37 


7 


4 


-1/10080 


80 


24 


3 


367/19353600 


123 


10 


8 


173/604800 


38 


8 


4 


17/40320 


81 


25 


3 


473/4838400 


124 


11 


8 


5/24192 


39 


6 


5 


-53/60480 


82 


26 


3 


41/215040 


125 


12 


8 


1/60480 


40 


7 


5 


-19/13440 


83 


27 


3 


211/1209600 


126 


13 


8 


61/403200 


41 


8 


5 


-1/5040 


84 


28 


3 


89/1209600 


127 


14 


8 


-1/151200 


42 


24 


1 





85 


29 


3 


1/100800 










43 


24 


2 





86 


30 


3 


79/967680 











26 



References 

[1] M. Abramowitz and LA. Stegun. Handbook of Mathematical Functions. 
Dover, 1965. 

[2] H.F. Baker. Alternant and continuous groups. Proc. London Math. Soc. 
(Second series), 3:24-47, 1905. 

[3] S. Blanes and F. Casas. On the convergence and optimization of the Baker- 
Campbell-Hausdorff formula. Linear Algebra AppL, 378:135-158, 2004. 

[4] S. Blanes, F. Casas, J. A. Oteo, and J. Ros. Magnus and Fer expansions for 
matrix differential equations: the convergence problem. J. Phys. A: Math. 
Gen., 22:259-268, 1998. 

[5] S. Blanes, F. Casas, J. A. Oteo, and J. Ros. Magnus expansion: mathe- 
matical study and physical applications, 2008. Work in progress. 

[6] A. Bose. Dynkin's method of computing the terms of the Baker-Campbell- 
Hausdorff series. J. Math. Phys., 30:2035-2037, 1989. 

[7] N. Bourbaki. Lie groups and Lie algebras. Chapters 1-3, Springer, 1989. 

[8] J. E. Campbell. On a law of combination of operators. Proc. London Math. 
Soc, 29:14-32, 1898. 

[9] F. Casas. Sufficient conditions for the convergence of the Magnus expan- 
sion. J. Phys. A: Math. Theor., 40:15001-15017, 2007. 

[10] A. J. Dragt and J. M. Finn. Lie series and invariant functions for analytic 
symplectic maps. J. Math. Phys., 17:2215-2227, 1976. 

[11] A. Diir. Mobius Functions, Incidence Algebras and Power-Series Repre- 
sentations, volume 1202 of Lecture Notes in Mathematics. Springer- Verlag, 
Berlin, 1986. 

[12] E.B. Dynkin. Evaluation of the coefficients of the Campbell-Hausdorff 
formula. Dokl. Akad. Nauk. SSSR, 57:323-326, 1947. 

[13] K. Goldberg. The formal power series for log(e^'e^). Duke Math. J., 23:13- 
21, 1956. 

[14] V.V. Gorbatsevich, A.L. Onishchik, and E.B. Vinberg. Foundations of Lie 
Theory and Lie Transformation Groups. Springer, 1997. 

[15] R. Grossman and R.G. Larson. Hopf-algebraic structure of families of trees. 
J. Algebra, 126:184-210, 1989. 

[16] E. Hairer, Ch. Lubich, and G. Wanner. Geometric Numerical Integra- 
tion. Structure-Preserving Algorithms for Ordinary Differential Equations. 
Springer- Verlag, Second edition, 2006. 



27 



[17] F. Hausdorff. Die symbolische Exponential Formel in der Gruppen theorie. 
Leipziger Ber., 58:19-48, 1906. 

[18] M.E. Hoffman. Combinatorics of rooted trees and Hopf algebras. Trans. 
Amer. Math. Soc, 355:3795-3811, 2003. 

[19] A. Iserles, H. Z. Munthe-Kaas, S. P. N0rsett, and A. Zanna. Lie-group 
methods. Acta Numerica, 9:215-365, 2000. 

[20] A. Iserles and S. P. N0rsett. On the solution of linear differential equations 
in Lie groups. Phil. Trans. Royal Soc. A, 357:983-1019, 1999. 

[21] N. Jacobson. Lie Algebras. Dover, 1979. 

[22] M. Kolsrud. Maximal reductions in the Baker-Hausdorff formula. J. Math. 
Phys., 34(l):270-285, 1993. 

[23] P.-V. Koseleff. Calcul formel pour les methodes de Lie en mecanique hamil- 
tonienne. PhD thesis, Ecole Polytechnique, 1993. 

[24] K. Kumar. On expanding the exponential. J. Math. Phys., 6:1928-1934, 
1965. 

[25] M. Lothaire. Combinatorics on Words, volume 17 of Encyclopedia of Math- 
ematics. Addison- Wesley, Reading, 1983. 

[26] W. Magnus. On the exponential solution of differential equations for a 
linear operator. Comm. Pure and Appl. Math., VII:649-673, 1954. 

[27] R.I. McLachlan and R. Quispel. Splitting methods. Acta Numerica, 
11:341-434, 2002. 

[28] J. Michel. Bases des algebres de Lie et serie de Hausdorff. Seminaire 
Dubreil. Algebre, 27(6): 1-9, 1974. 

[29] P.C. Moan. On backward error analysis and Nekhoroshev stability in the 
numerical analysis of conservative systems of DDEs. PhD thesis. Univer- 
sity of Cambridge, 2002. 

[30] A. Murua. The Hopf algebra of rooted trees, free Lie algebras, and Lie 
series. Found. Comp. Math., 6:387-426, 2006. 

[31] M. Newman, W. So, and R.C. Thompson. Convergence domains for the 
Campbell-Baker-Hausdorff formula. Linear Multilinear Algebra, 24:301- 
310, 1989. 

[32] M. Newman and R.C. Thompson. Numerical values of Goldberg's coeffi- 
cients in the series for log(e'"e2'). Math. Comp., 48:265-271, 1987. 

[33] J. A. Oteo. The Baker-Campbell-Hausdorff formula and nested commutator 
identities. J. Math. Phys., 32(2):419-424, 1991. 



28 



[34] M. Postnikov. Lie Groups and Lie Algebras. Semester V of Lectures in 
Geometry. URSS Publishers, 1994. 

[35] M.W. Reinsch. A simple expression for the terms in the Baker-Campbell- 
Hausdorff series. J. Math. Phys., 41:2434-2442, 2000. 

[36] C. Reutenauer. Free Lie algebras, volume 7 of London Math. Soc. mono- 
graphs (new series). Oxford University Press, 1993. 

[37] R.D. Richtmyer and S. Greenspan. Expansion of the Campbell-Baker- 
Hausdorff formula by computer. Gomm. Pure and Appl. Math., 18:107- 
108, 1965. 

[38] F. Rouviere. Espaces symetriques et methode de Kashiwara-Vergne. Ann. 
Sci. Ecole Norm. Sup., 19:553-581, 1986. 

[39] J. M. Sanz-Serna and M. P. Calvo. Numerical Hamiltonian Problems. 
AMMC 7. Chapman & Hall, 1994. 

[40] A.T. Sornborger and E.D. Stewart. Higher-order methods for simulations 
on quantum computers. Phys. Rev. A, 60:1956-1965, 1999. 

[41] M. Suzuki. On the convergence of exponential operators — the Zassenhaus 
formula, BCH formula and systematic approximants. Commun. Math. 
Phys., 57:193-200, 1977. 

[42] R.C. Thompson. Cyclic relations and the Goldberg coefficients in the 
Campbell-Baker-Hausdorff formula. Proc. Amer. Math. Soc, 86:12-14, 
1982. 

[43] R.C. Thompson. Proof of a conjectured exponential formula. Linear Mul- 
tilinear Algebra, 19:187-197, 1986. 

[44] R.C. Thompson. Special cases of a matrix exponential formula. Linear 
Algebra Appl, 107:283-292, 1988. 

[45] R.C. Thompson. Convergence proof for Goldberg's exponential series. Lin- 
ear Algebra Appl, 121:3-7, 1989. 

[46] M. Torres- Torriti and H. Michalska. A software package for Lie algebraic 
computations. SIAM Review, 47:722-745, 2005. 

[47] V. S. Varadarajan. Lie Groups, Lie Algebras, and Their Representations. 
GTM 102. Springer- Verlag, 1984. 

[48] X.G. Viennot. Algebres de Lie Libres et monoides libres, volume 691 of 
Lecture Notes in Mathematics. Springer, Berlin, 1978. 

[49] G.H. Weiss and A. A. Maradudin. The Baker-Hausdorff formula and a 
problem in Crystal Physics. J. Math. Phys., 3:771-777, 1962. 

[50] R.M. Wilcox. Exponential operators and parameter differentiation in quan- 
tum physics. J. Math. Phys., 8:962-982, 1967. 



29 



[51] V.A. Yakubovich and V.M. Starzhinskii. Linear Differential Equations 
with Periodic Coefficients. John Wiley & Sons, 1975. 

[52] H. Yoshida. Construction of higher order symplectic integrators. Phys. 
Lett. A, 150:262-268, 1990. 



30 



